Goal: Test common assumptions about firefly size (elytron length and body mass) over the season and across locations in P. pyralis fireflies

Study design:

Step 1: Ready the workspace

#clear all inputs
rm(list = ls())

#Check for necessary packages
#list.of.packages <- c("ggplot2", 
#                      "tidyr",
#                      "dplyr",
#                      "Hmisc",
#                      "qqplotr",
#                      "ggthemes",
#                      "fabricatr",
#                      "gridExtra",
#                      "grid",
#                      "kableExtra",
#                      "sjPlot",
#                      "cowplot",
#                      "ggpubr",
#                      "patchwork",
#                      "magick", #so you can use the savekable function
#                      "webshot", #so you can use the savekable function
#                      "lme4",
#                      "RcppEigen")

#should install packages that you don't have
#new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
#if(length(new.packages)) install.packages(new.packages)

#Load the package that contains the Cox Propotional Hazard Function
library(ggplot2)
library(tidyr)
library(dplyr)
library(Hmisc)
library(qqplotr)
library(ggthemes)
library(fabricatr)
library(gridExtra)
library(grid)
library(kableExtra)
library(sjPlot)
library(cowplot)
library(ggpubr)
library(patchwork)

Step 2: Import the data

#This imports the data from a CSV file and creates the data frame PrTC
Adult_Final<-read.table("2021_Adult_Firefly_Survival_binned3.csv", header=TRUE, sep=",", dec=".",na.strings=".") 

#This sets up the color palette using tableau20: https://jrnold.github.io/ggthemes/reference/tableau_color_pal.html
site_colors <- c("USA: Montour Co, Bucknell Natural Area" = "#4E79A7", 
                 "USA: Union Co, Bucknell Farm" = "#F28E2B", 
                 "USA: Union Co, Bucknell Ropes Course" = "#E15759")

Data Used:

  • 838 fireflies with elytra length (calipers), body mass (analytical balance), and survival data.

Abbreviations used:

  • BNA = Bucknell Natural Area
  • BUF = Bucknell University Farm
  • FDBCC = Forrest D. Brown Conference Center

Step 3: Initial Analysis

Elytron length

Q1. What does the length distribution look like overall?

  • Raw data (left) vs log-transformed (right)
#generate distribution of raw data
z <- ggplot(Adult_Final, aes(x=ElytralLength)) +
  geom_histogram() +
  xlab("Elytron length (mm)") +
  ylab("N fireflies") +
  theme(plot.title = element_text(hjust = 0.5))

#generate log distribution plot
y <- ggplot(Adult_Final, aes(x=log(ElytralLength))) +
  geom_histogram() +
  xlab("log(Elytron length) (mm)") +
  ylab("N fireflies") +
  theme(plot.title = element_text(hjust = 0.5))

grid.arrange(z, y, nrow= 1)

Q2. …by location?

  • Data summary (untransformed):
#generate tibble
t<- Adult_Final %>% 
  group_by(Location) %>% 
  summarise(N=n(), mean_length=round(mean(ElytralLength),2), SD = round(sd(ElytralLength), 2))

#generate nice looking table
kbl(t, col.names = c("Location", "n", "Mean elytron length (mm)", "SD"), align = "lccc") %>%
  kable_classic(full_width = F, html_font = "Cambria")
Location n Mean elytron length (mm) SD
USA: Montour Co, Bucknell Natural Area 298 9.11 0.96
USA: Union Co, Bucknell Farm 70 8.90 1.08
USA: Union Co, Bucknell Ropes Course 470 8.94 0.95
  • Raw data, three different ways:
#generate labels for sites that include sample size
site_labels = c(paste("BNA", " (n = ", t$N[1], ")", sep=""), 
                paste0("BUF", " (n = ", t$N[2], ")"), 
                paste0("FDBCC",  " (n = ", t$N[3], ")"))

#generate text placeholder
a <- ggplot(Adult_Final, aes(x=ElytralLength, fill=Location)) +
  geom_histogram(alpha= 0.5, position="identity") +
  scale_fill_manual(name = "Collection Location", labels = site_labels, values = site_colors)
  
legend.a <- cowplot::get_legend(a)

#generate histogram
b <- ggplot(Adult_Final, aes(x=ElytralLength, fill=Location)) +
  geom_histogram(alpha= 0.5, position="identity") +
  xlab("Elytron length (mm)") +
  ylab("N fireflies") +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none") +
  scale_fill_manual(name = "Location", labels = site_labels, values = site_colors)

#generate violin plot
c <- ggplot(Adult_Final, aes(x=Location, y=ElytralLength, fill=Location)) +
  geom_violin(alpha = 0.9) +
  theme_classic() +
  ylab("Elytron length (mm)") +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none", axis.text.x=element_blank(), axis.title.x=element_blank()) +
  geom_boxplot(width=0.2, alpha= 0.5, fill = "white") +
  scale_fill_manual(values = site_colors)

#generate boxplot
d <- ggplot(Adult_Final, aes(x=Location, y=ElytralLength, fill=Location)) +
  geom_boxplot(alpha = 0.9) +
  geom_jitter(size = 0.1)+
  theme_classic() +
  ylab("Elytron length (mm)") +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none", axis.text.x=element_blank(), axis.title.x=element_blank()) +
  scale_fill_manual(values = site_colors)

grid.arrange(legend.a,b,c,d, nrow=2)

Q3. …by seasonday

  • Data summary (untransformed):
#generate tibble
t <- Adult_Final %>%
  group_by(Location, SeasonDay) %>%
  summarise(N=n(), mean_length=round(mean(ElytralLength),2), SD = round(sd(ElytralLength),2)) %>%
  arrange(SeasonDay)

#make table nice in kable
kbl(t, col.names = c("Location", "Season Day", "n", "Mean elytron length (mm)", "SD"), align = "lcccc") %>%
  kable_classic(full_width = F, html_font = "Cambria")
Location Season Day n Mean elytron length (mm) SD
USA: Montour Co, Bucknell Natural Area 9 92 9.55 0.91
USA: Union Co, Bucknell Farm 13 47 8.76 1.09
USA: Montour Co, Bucknell Natural Area 15 120 8.71 0.98
USA: Union Co, Bucknell Ropes Course 17 231 8.96 0.94
USA: Montour Co, Bucknell Natural Area 22 86 9.21 0.73
USA: Union Co, Bucknell Ropes Course 23 95 8.93 1.01
USA: Union Co, Bucknell Farm 30 23 9.21 1.01
USA: Union Co, Bucknell Ropes Course 31 90 8.94 0.93
USA: Union Co, Bucknell Ropes Course 36 54 8.90 0.97
  • Raw data, three different ways:
#generate labels for seasondays that include sample size
t$SeasonDayLabel = paste(t$SeasonDay, " (", t$N, ")", sep ="")
  
seasonday_label <- c(
  '9' = t$SeasonDayLabel[1],
  '13' = t$SeasonDayLabel[2],
  '15' = t$SeasonDayLabel[3],
  '17' = t$SeasonDayLabel[4],
  '22' = t$SeasonDayLabel[5],
  '23' = t$SeasonDayLabel[6],
  '30' = t$SeasonDayLabel[7],
  '31' = t$SeasonDayLabel[8],
  '36' = t$SeasonDayLabel[9]
  )

#generate legend
e <- ggplot(Adult_Final, aes(x=ElytralLength, fill=Location)) +
  geom_histogram() +
  xlab("Elytron length (mm)") +
  ylab("N fireflies") +
  theme(plot.title = element_text(hjust = 0.5)) +
  facet_wrap(~SeasonDay, labeller = as_labeller(seasonday_label)) +
  scale_fill_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759"))

legend.e <- cowplot::get_legend(e)

#generate histogram
f <- ggplot(Adult_Final, aes(x=log(ElytralLength), fill=Location)) +
  geom_histogram() +
  xlab("Elytron length (mm)") +
  ylab("N fireflies") +
  ggtitle("SeasonDay (n)") +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none") +
  facet_wrap(~SeasonDay, labeller = as_labeller(seasonday_label), ncol=5) +
  scale_fill_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759"))

#generate violin plot by seasonday
g <- ggplot(Adult_Final, aes(x=as.factor(SeasonDay), y=ElytralLength, fill = Location)) +
  geom_violin(alpha = 0.9) +
  theme_classic() +
  ylab("Elytron length (mm)") +
  xlab("Season Day") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90), legend.position = "none") +
  geom_boxplot(width=0.2, fill = "white") +
  scale_fill_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759"))

#generate violin plot by site
h <- ggplot(Adult_Final, aes(x=as.factor(SeasonDay), y=ElytralLength, fill = Location)) +
  geom_violin(alpha = 0.9) +
  theme_classic() +
  ylab("Elytron length (mm)") +
  xlab("Season Day") +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none", axis.text.x = element_text(angle = 90)) +
  facet_wrap(~Location, labeller = labeller(Location = 
    c("USA: Montour Co, Bucknell Natural Area" = "BNA",
      "USA: Union Co, Bucknell Farm" = "BUF",
      "USA: Union Co, Bucknell Ropes Course" = "FDBCC"))) +
  geom_boxplot(width=0.2, fill = "white") +
  scale_fill_manual(values = c("#4E79A7", "#F28E2B", "#E15759"))

lay1 <- rbind(c(2,2,2,1))
lay2 <- rbind(c(1,1,2,2,2))

g1 <- arrangeGrob(legend.e, f, layout_matrix = lay1)
g2 <- arrangeGrob(g,h, layout_matrix = lay2)
 
grid.arrange(g1, g2, nrow=2)

Q4. Are the distributions within each Location * SeasonDay normal?

Testing the raw data:

  • Q-Q plot for the raw data (left) and Shapiro test (bottom right)
#qq plots to assess normality
q <- ggplot(Adult_Final, aes(sample=ElytralLength, color=Location)) +
  stat_qq_band(bandType = "pointwise", fill = "#8DA0CB", alpha = 0.4) + 
  stat_qq_line(color = "#8DA0CB") +
  stat_qq_point() +
  ggtitle("Season Day") +
  facet_wrap(~SeasonDay) +
  scale_color_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759")) +
  theme(legend.position = "none")

#statistically test
shapiro.table <- NULL
for (day in levels(factor(Adult_Final$SeasonDay))){
  #subset table
  table.sub <- Adult_Final[Adult_Final$SeasonDay == day,]
  test.result <- (shapiro.test(table.sub$ElytralLength))
  info <- data.frame(SeasonDay = day, W = round(test.result[1]$statistic, 2), p = round(test.result[2]$p.value, 2))
  shapiro.table <- rbind(shapiro.table, info)
}

#create kable filename
kable_filename = paste("Elytron_length_shapiro", ".png", sep="")

#generate nice table with kable
kbl(shapiro.table, row.names = FALSE, col.names = c("Season Day", "W", "P"), align = "ccc")  %>%
  kable_classic(full_width = F, html_font = "Cambria") %>% 
  save_kable(file = kable_filename, bs_theme = "cerulean", zoom = 7) 

#make big plot with table
#read in kable table as image
r <- ggdraw() +
  draw_image(kable_filename)

#stack the legend on top of the table
g3 <- arrangeGrob(legend.e, r, nrow =2)

#make the column width layout
shapiro_layout <- rbind(c(1,1,1,2))

#plot the graph next to the legend and table
grid.arrange(q, g3, nrow = 1,  layout_matrix = shapiro_layout)

Conclusions:

  • p >> 0.05 for elytron length distributions from each collection day, and thus, cannot reject null (normality).

This suggests that log transformation is not necessary to satisfy statistical assumptions of a linear relationship between the response variable (elytron length) and explanatory variables (season day, location) when looking at how these variables are related. However, log transformation is the standard when comparing body size metrics due to allometric scaling expectations. Thus, we will proceed with testing whether log transformation of elytron length alters expected linear relationships among variables (specifically body mass and survival). The fits may be no different or even better.

Q5. Does log-transformation change things?

Testing log-transformed data:

  • Q-Q plot for the log-transformed data (left) and Shapiro test (bottom right)
#qq plots to assess normality
q <- ggplot(Adult_Final, aes(sample=log(ElytralLength), color=Location)) +
  stat_qq_band(bandType = "pointwise", fill = "#8DA0CB", alpha = 0.4) + 
  stat_qq_line(color = "#8DA0CB") +
  stat_qq_point() +
  ggtitle("Season Day") +
  facet_wrap(~SeasonDay) +
  scale_color_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759")) +
  theme(legend.position = "none")

#statistically test
shapiro.table <- NULL
for (day in levels(factor(Adult_Final$SeasonDay))){
  #subset table
  table.sub <- Adult_Final[Adult_Final$SeasonDay == day,]
  test.result <- (shapiro.test(log(table.sub$ElytralLength)))
  info <- data.frame(SeasonDay = day, W = round(test.result[1]$statistic, 2), p = round(test.result[2]$p.value, 2))
  shapiro.table <- rbind(shapiro.table, info)
}

#create kable filename
kable_filename = paste("log_Elytron_length_shapiro", ".png", sep="")

#generate nice table with kable
kbl(shapiro.table, row.names = FALSE, col.names = c("Season Day", "W", "P"), align = "ccc")  %>%
  kable_classic(full_width = F, html_font = "Cambria") %>% 
  save_kable(file = kable_filename, bs_theme = "cerulean", zoom = 7) 

#make big plot with table
#read in kable table as image
r <- ggdraw() +
  draw_image(kable_filename)

#stack the legend on top of the table
g3 <- arrangeGrob(legend.e, r, nrow =2)

#make the column width layout
shapiro_layout <- rbind(c(1,1,1,2))

#plot the graph next to the legend and table
grid.arrange(q, g3, nrow = 1,  layout_matrix = shapiro_layout)

Conclusions:

  • p < 0.05 for log elytron length distributions from one collection day, and thus, reject normality. Two others (day 13 and day 23) are close (p = 0.08 and p = 0.06, respectively). Day 13 likely has smaller sample size (BUF). Thus, most days are consistent with normal distribution.
  • Since the distributions more frequently depart from normality, we will use untransformed (raw) elytron length in analyses investigating relationships with Season Day and Location, but will also investigate the effects of log transformation in elytron length by body mass comparisons. Perhaps transformation will not make a big difference. Most sample sizes are > 30, meaning they are large enough for parametric tests anyway.

Q6. Do seasonday and location predict elytron length?

  • Summary of effects:
Elytral<-lm(ElytralLength~SeasonDay*Location, data=Adult_Final)
summary(Elytral)
## 
## Call:
## lm(formula = ElytralLength ~ SeasonDay * Location, data = Adult_Final)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.98747 -0.62992  0.00396  0.64271  2.91253 
## 
## Coefficients:
##                                                        Estimate Std. Error
## (Intercept)                                             9.48258    0.17763
## SeasonDay                                              -0.02430    0.01112
## LocationUSA: Union Co, Bucknell Farm                   -1.06859    0.34167
## LocationUSA: Union Co, Bucknell Ropes Course           -0.48694    0.23421
## SeasonDay:LocationUSA: Union Co, Bucknell Farm          0.05070    0.01821
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course  0.02206    0.01279
##                                                        t value Pr(>|t|)    
## (Intercept)                                             53.385  < 2e-16 ***
## SeasonDay                                               -2.186  0.02910 *  
## LocationUSA: Union Co, Bucknell Farm                    -3.128  0.00182 ** 
## LocationUSA: Union Co, Bucknell Ropes Course            -2.079  0.03792 *  
## SeasonDay:LocationUSA: Union Co, Bucknell Farm           2.784  0.00550 ** 
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course   1.724  0.08505 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9639 on 832 degrees of freedom
## Multiple R-squared:  0.01736,    Adjusted R-squared:  0.01145 
## F-statistic:  2.94 on 5 and 832 DF,  p-value: 0.01223
  • ANOVA(model) results:
kbl(anova(Elytral), digits = 2)  %>%
  kable_classic(full_width = F, html_font = "Cambria")
Df Sum Sq Mean Sq F value Pr(>F)
SeasonDay 1 2.66 2.66 2.87 0.09
Location 2 3.71 1.86 2.00 0.14
SeasonDay:Location 2 7.28 3.64 3.92 0.02
Residuals 832 773.03 0.93 NA NA
  • For Elytron length, only the interaction term (SeasonDay:Location) is significant (p = 0.02).

Conclusion

  • Use log-transformed elytron length in downstream analysis when using elytron length as an explanatory variable rather than a response variable with respect to variation in body mass due to allometric scaling.

  • Season Day doesn’t really matter for elytron length (this is expected if elytron length is fixed for the life of adult).

  • Location doesn’t really matter (caveats: sampling efforts at sites were not equal and non-randomly distributed across the season)

  • SeasonDay:Location is significant. This is apparent if you look at the violin plots for Season Day grouped by Location. There is a downward trend in BNA, but not in the others. However, the sample size (replicate sampling events at a particular site across the season) is small (BNA = sampled 3 times in the first half of the season, BUF = sampled twice, FDBCC = sampled 4 times, all in second half of season).

  • Side Note: Vencl and Carlson (1998) (https://www.proquest.com/docview/757065622?accountid=9784) found a bimodal distribution in elytron length at their sites on Long Island across season (July), no effect of year (3 field seasons) so these were all pooled in the distribution displayed in the manuscript. They log-transformed elytron length for their analysis, and the plot displays the untransformed data.

  • Plots:

    • Blue: BNA
    • Red: BUF
    • Gold: FDBCC
Overall plot
#plot with all locations over the season
p_elytral_length_by_season_day <- ggplot(Adult_Final, aes(x=as.numeric(SeasonDay), y=ElytralLength, color = Location)) +
  geom_point(size = 1, alpha = 0.7) +
  theme_classic() +
  ylab("Elytron length (mm)") +
  xlab("Season Day") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90), legend.position = "none") +
  scale_color_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759"))

p_elytral_length_by_season_day

Investigating individual locations
BNA
  • Stats:
#BNA
Elytral_BNA<-lm(ElytralLength~SeasonDay, data=Adult_Final[Adult_Final$Location == "USA: Montour Co, Bucknell Natural Area",])
summary(Elytral_BNA)
## 
## Call:
## lm(formula = ElytralLength ~ SeasonDay, data = Adult_Final[Adult_Final$Location == 
##     "USA: Montour Co, Bucknell Natural Area", ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.79806 -0.62992  0.07903  0.66613  2.41613 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   9.4826     0.1758  53.937   <2e-16 ***
## SeasonDay    -0.0243     0.0110  -2.209    0.028 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.954 on 296 degrees of freedom
## Multiple R-squared:  0.01621,    Adjusted R-squared:  0.01289 
## F-statistic: 4.878 on 1 and 296 DF,  p-value: 0.02797
p_elytral_length_by_season_day_BNA <- ggplot(Adult_Final[Adult_Final$Location == "USA: Montour Co, Bucknell Natural Area",], aes(x=as.numeric(SeasonDay), y=ElytralLength)) +
  geom_point(size = 1, alpha = 0.7, color = "#4E79A7") +
  theme_classic() +
  ylab("Elytron length (mm)") +
  xlab("Season Day") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90), legend.position = "none") +
  geom_smooth(method = "lm", se=FALSE, color = "gray23") +
  stat_regline_equation(label.x = 17, label.y = 11.4, aes(label = ..eq.label..)) +
  stat_regline_equation(label.x = 17, label.y = 11.0, aes(label = ..rr.label..)) +
  ggtitle("BNA")

#p_elytral_length_by_season_day_BNA
BUF
  • Stats:
Elytral_BUF<-lm(ElytralLength~SeasonDay, data=Adult_Final[Adult_Final$Location == "USA: Union Co, Bucknell Farm",])
summary(Elytral_BUF)
## 
## Call:
## lm(formula = ElytralLength ~ SeasonDay, data = Adult_Final[Adult_Final$Location == 
##     "USA: Union Co, Bucknell Farm", ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.27723 -0.61695  0.02277  0.75613  2.02277 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.41399    0.32250  26.090   <2e-16 ***
## SeasonDay    0.02640    0.01594   1.656    0.102    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.065 on 68 degrees of freedom
## Multiple R-squared:  0.03877,    Adjusted R-squared:  0.02463 
## F-statistic: 2.743 on 1 and 68 DF,  p-value: 0.1023
p_elytral_length_by_season_day_BUF <- ggplot(Adult_Final[Adult_Final$Location == "USA: Union Co, Bucknell Farm",], aes(x=as.numeric(SeasonDay), y=ElytralLength)) +
  geom_point(size = 1, alpha = 0.7, color = "#F28E2B") +
  theme_classic() +
  ylab("Elytron length (mm)") +
  xlab("Season Day") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90), legend.position = "none") +
  ggtitle("BUF")

#p_log_elytral_length_by_season_day_BUF
FDBCC
  • Stats:
Elytral_FDBCC<-lm(ElytralLength~SeasonDay, data=Adult_Final[Adult_Final$Location == "USA: Union Co, Bucknell Ropes Course",])
summary(Elytral_FDBCC)
## 
## Call:
## lm(formula = ElytralLength ~ SeasonDay, data = Adult_Final[Adult_Final$Location == 
##     "USA: Union Co, Bucknell Ropes Course", ])
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.98747 -0.63048 -0.03074  0.62789  2.91253 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.995642   0.151183  59.502   <2e-16 ***
## SeasonDay   -0.002245   0.006267  -0.358     0.72    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9546 on 468 degrees of freedom
## Multiple R-squared:  0.0002742,  Adjusted R-squared:  -0.001862 
## F-statistic: 0.1283 on 1 and 468 DF,  p-value: 0.7203
p_elytral_length_by_season_day_FDBCC <- ggplot(Adult_Final[Adult_Final$Location == "USA: Union Co, Bucknell Ropes Course",], aes(x=as.numeric(SeasonDay), y=ElytralLength)) +
  geom_point(size = 1, alpha = 0.7, color = "#E15759") +
  theme_classic() +
  ylab("Elytron length (mm)") +
  xlab("Season Day") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90), legend.position = "none")  +
  ggtitle("FDBCC")

#p_log_elytral_length_by_season_day_FDBCC
Plots
p_elytral_length_by_season_day_BNA + p_elytral_length_by_season_day_BUF + p_elytral_length_by_season_day_FDBCC

Conclusions:

  • elytron length decreases over the season at BNA (p = 0.03), with 3 sampling days
  • doesn’t have an effect at BUF (p = 0.1), but only 2 sampling days
  • doesn’t have an effect at FDBCC (p = 0.7), with 4 sampling days
  • additional replicated sampling days at each site, preferably distributed across the season, are required to further investigate this trend

Body mass

Q1. What does the mass distribution look like overall?

  • Raw data (left) vs log-transformed (right)
#generate distribution of raw data
z <- ggplot(Adult_Final, aes(x=Mass)) +
  geom_histogram() +
  xlab("Body mass (g)") +
  ylab("N fireflies") +
  theme(plot.title = element_text(hjust = 0.5))

#generate log distribution plot
y <- ggplot(Adult_Final, aes(x=log(Mass))) +
  geom_histogram() +
  xlab("log(Body mass) (g)") +
  ylab("N fireflies") +
  theme(plot.title = element_text(hjust = 0.5))

grid.arrange(z, y, nrow= 1)

Q2. …by location?

  • Data summary (untransformed):
#generate tibble 
t<- Adult_Final %>% 
  group_by(Location) %>% 
  summarise(N=n(), mean_mass=round(mean(Mass),4), SD = round(sd(Mass), 4))

#generate nice looking table
kbl(t, col.names = c("Location", "n", "Mean body mass (g)", "SD"), align = "lccc") %>%
  kable_classic(full_width = F, html_font = "Cambria")
Location n Mean body mass (g) SD
USA: Montour Co, Bucknell Natural Area 298 0.0323 0.0092
USA: Union Co, Bucknell Farm 70 0.0326 0.0088
USA: Union Co, Bucknell Ropes Course 470 0.0312 0.0083
  • Log-transformed data, three different ways:
#generate sample size labels
site_labels = c(paste("BNA", " (N=", t$N[1], ")", sep=""), 
                paste0("BUF", " (N=", t$N[2], ")"), 
                paste0("FDBCC",  " (N=", t$N[3], ")"))

#generate text placeholder
i <- ggplot(Adult_Final, aes(x=log(Mass), fill=Location)) +
  geom_histogram(alpha= 0.5, position="identity") +
  scale_fill_manual(name = "Collection Location", labels = site_labels, values = site_colors)
  
legend.i <- cowplot::get_legend(i)

#generate histogram with log mass
j <- ggplot(Adult_Final, aes(x=log(Mass), fill=Location)) +
  geom_histogram(alpha= 0.5, position="identity") +
  xlab("log(Body mass) (g)") +
  ylab("N fireflies") +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none") +
  scale_fill_manual(name = "Location", labels = site_labels, values = site_colors)

#generate violin plot with log mass
k <- ggplot(Adult_Final, aes(x=Location, y=log(Mass), fill=Location)) +
  geom_violin(alpha = 0.9) +
  theme_classic() +
  ylab("log(Body mass) (g)") +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none", axis.text.x=element_blank(), axis.title.x=element_blank()) +
  geom_boxplot(width=0.2, alpha= 0.5, fill = "white") +
  scale_fill_manual(values = site_colors)

#generate boxplot
l <- ggplot(Adult_Final, aes(x=Location, y=log(Mass), fill=Location)) +
  geom_boxplot(alpha = 0.9) +
  geom_jitter(size = 0.1)+
  theme_classic() +
  ylab("log(Body mass) (g)") +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none", axis.text.x=element_blank(), axis.title.x=element_blank()) +
  scale_fill_manual(values = site_colors)

grid.arrange(legend.i,j,k,l, nrow=2)

Q3. …by seasonday

  • Data summary (untransformed):
#generate tibble
t <- Adult_Final %>%
  group_by(Location, SeasonDay) %>%
  summarise(N=n(), mean_mass=round(mean(Mass),4), SD = round(sd(Mass),4)) %>%
  arrange(SeasonDay)

#make table nice in kable
kbl(t, col.names = c("Location", "Season Day", "n", "Mean body mass (g)", "SD"), align = "lcccc") %>%
  kable_classic(full_width = F, html_font = "Cambria")
Location Season Day n Mean body mass (g) SD
USA: Montour Co, Bucknell Natural Area 9 92 0.0337 0.0092
USA: Union Co, Bucknell Farm 13 47 0.0325 0.0095
USA: Montour Co, Bucknell Natural Area 15 120 0.0313 0.0099
USA: Union Co, Bucknell Ropes Course 17 231 0.0321 0.0087
USA: Montour Co, Bucknell Natural Area 22 86 0.0320 0.0081
USA: Union Co, Bucknell Ropes Course 23 95 0.0312 0.0088
USA: Union Co, Bucknell Farm 30 23 0.0330 0.0074
USA: Union Co, Bucknell Ropes Course 31 90 0.0292 0.0071
USA: Union Co, Bucknell Ropes Course 36 54 0.0306 0.0072
  • Log-transformed data, three different ways:
#generate labels for seasondays that include sample size
t$SeasonDayLabel = paste(t$SeasonDay, " (", t$N, ")", sep ="")

seasonday_label <- c(
  '9' = t$SeasonDayLabel[1],
  '13' = t$SeasonDayLabel[2],
  '15' = t$SeasonDayLabel[3],
  '17' = t$SeasonDayLabel[4],
  '22' = t$SeasonDayLabel[5],
  '23' = t$SeasonDayLabel[6],
  '30' = t$SeasonDayLabel[7],
  '31' = t$SeasonDayLabel[8],
  '36' = t$SeasonDayLabel[9]
  )

#generate histogram
m <- ggplot(Adult_Final, aes(x=log(Mass), fill=Location)) +
  geom_histogram() +
  xlab("log(Body mass) (g)") +
  ylab("N fireflies") +
  theme(plot.title = element_text(hjust = 0.5)) +
  facet_wrap(~SeasonDay, labeller = as_labeller(seasonday_label)) +
  scale_fill_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759")) 

#generate legend
legend.m <- cowplot::get_legend(m)

#generate histogram
n <- ggplot(Adult_Final, aes(x=log(Mass), fill=Location)) +
  geom_histogram() +
  xlab("log(Body mass) (g)") +
  ylab("N fireflies") +
  ggtitle("Season Day (n)") +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none") +
  facet_wrap(~SeasonDay, labeller = as_labeller(seasonday_label), ncol = 5) +
  scale_fill_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759")) 

#generate violin plot by seasonday
o <- ggplot(Adult_Final, aes(x=as.factor(SeasonDay), y=log(Mass), fill = Location)) +
  geom_violin(alpha = 0.9) +
  theme_classic() +
  ylab("log(Body mass) (g)") +
  xlab("Season Day") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90), legend.position = "none") +
  geom_boxplot(width=0.2, fill = "white") +
  scale_fill_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759"))

#generate violin plot by site
p <- ggplot(Adult_Final, aes(x=as.factor(SeasonDay), y=log(Mass), fill = Location)) +
  geom_violin(alpha = 0.9) +
  theme_classic() +
  ylab("log(Body mass) (g)") +
  xlab("Season Day") +
  theme(plot.title = element_text(hjust = 0.5), legend.position = "none", axis.text.x = element_text(angle = 90)) +
  facet_wrap(~Location, labeller = labeller(Location = 
    c("USA: Montour Co, Bucknell Natural Area" = "BNA",
      "USA: Union Co, Bucknell Farm" = "BUF",
      "USA: Union Co, Bucknell Ropes Course" = "FDBCC"))) +
  geom_boxplot(width=0.2, fill = "white") +
  scale_fill_manual(values = c("#4E79A7", "#F28E2B", "#E15759"))

lay1 <- rbind(c(2,2,2,1))
lay2 <- rbind(c(1,1,2,2,2))

g3 <- arrangeGrob(legend.m, n, layout_matrix = lay1)
g4 <- arrangeGrob(o,p, layout_matrix = lay2)
 
grid.arrange(g3, g4, nrow=2)

Q4. Are the distributions within each Site * SeasonDay normal?

Testing the raw data:

  • Q-Q plots for the raw data (left) and Shapiro test (bottom right)
#qq plots to assess normality with raw data
q <- ggplot(Adult_Final, aes(sample=Mass, color=Location)) +
  stat_qq_band(bandType = "pointwise", fill = "#8DA0CB", alpha = 0.4) + 
  stat_qq_line(color = "#8DA0CB") +
  stat_qq_point() +
  facet_wrap(~SeasonDay) +
  scale_color_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759")) +
  theme(legend.position = "none")

#statistically test
shapiro.table <- NULL
for (day in levels(factor(Adult_Final$SeasonDay))){
  #subset table
  table.sub <- Adult_Final[Adult_Final$SeasonDay == day,]
  test.result <- (shapiro.test(table.sub$Mass))
  info <- data.frame(SeasonDay = day, W = round(test.result[1]$statistic, 2), p = round(test.result[2]$p.value, 2))
  shapiro.table <- rbind(shapiro.table, info)
}

#generate kable filename
kable_filename = paste("Body_mass_raw_shapiro", ".png", sep="")

#generate nice table with kable
kbl(shapiro.table, row.names = FALSE, col.names = c("Season Day", "W", "P"), align = "ccc")  %>%
  kable_classic(full_width = F, html_font = "Cambria") %>%
  save_kable(file = kable_filename, bs_theme = "cerulean", zoom = 7)

#read in kable table as image
r <- ggdraw() +
  draw_image(kable_filename)

#stack the legend on top of the table
g3 <- arrangeGrob(legend.e, r, nrow =2)

#make the column width layout
shapiro_layout <- rbind(c(1,1,1,2))

#plot the graph next to the legend and table
grid.arrange(q, g3, nrow = 1,  layout_matrix = shapiro_layout)

Conclusions: Not normal (5/9 days):

  • Day 13
  • Day 15
  • Day 17
  • Day 23
  • Day 31

However, sample size is large on these days - only Day 30 (not included in list of collections that violate normality) has sample size < 30. Next, test log-transformed data to see if improves.

Q5. Does log-transformation help?

Note: Log-transformation of body mass is necessary for the assumption of linearity in survival models (see Methods in main text)

Testing log-transformed data:

  • Q-Q plots for the log-transformed data (left) and Shapiro test (bottom right)
#qq plots to assess normality
q <- ggplot(Adult_Final, aes(sample=log(Mass), color=Location)) +
  stat_qq_band(bandType = "pointwise", fill = "#8DA0CB", alpha = 0.4) + 
  stat_qq_line(color = "#8DA0CB") +
  stat_qq_point() +
  facet_wrap(~SeasonDay) +
  scale_color_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759")) +
  theme(legend.position = "none")

#statistically test
shapiro.table <- NULL
for (day in levels(factor(Adult_Final$SeasonDay))){
  #subset table
  table.sub <- Adult_Final[Adult_Final$SeasonDay == day,]
  test.result <- (shapiro.test(log(table.sub$Mass)))
  info <- data.frame(SeasonDay = day, W = round(test.result[1]$statistic, 2), p = round(test.result[2]$p.value, 2))
  shapiro.table <- rbind(shapiro.table, info)
}

#generate kable table filename
kable_filename = paste("Body_mass_log_shapiro", ".png", sep="")
  
#generate nice table with kable
kbl(shapiro.table, row.names = FALSE, col.names = c("Season Day", "W", "P"), align = "ccc")  %>%
  kable_classic(full_width = F, html_font = "Cambria") %>%
  save_kable(file = kable_filename, bs_theme = "cerulean", zoom = 7)

#read in kable table as image
r <- ggdraw() +
  draw_image(kable_filename)

#stack the legend on top of the table
g3 <- arrangeGrob(legend.e, r, nrow =2)

#make the column width layout
shapiro_layout <- rbind(c(1,1,1,2))

#plot the graph next to the legend and table
grid.arrange(q, g3, nrow = 1,  layout_matrix = shapiro_layout)

Conclusions:

  • Yes - P >> 0.05 for log body mass distributions from each collection day, fail to reject normality.

Q6. Do seasonday and location predict body mass?

  • Summary of effects:
#no Elytral length; + interaction
Mass<-lm(log(Mass)~SeasonDay*Location, data=Adult_Final)
summary(Mass)
## 
## Call:
## lm(formula = log(Mass) ~ SeasonDay * Location, data = Adult_Final)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.74613 -0.18636  0.00805  0.18811  0.72058 
## 
## Coefficients:
##                                                          Estimate Std. Error
## (Intercept)                                            -3.4227082  0.0501670
## SeasonDay                                              -0.0034295  0.0031398
## LocationUSA: Union Co, Bucknell Farm                   -0.0684593  0.0964965
## LocationUSA: Union Co, Bucknell Ropes Course            0.0129793  0.0661484
## SeasonDay:LocationUSA: Union Co, Bucknell Farm          0.0052895  0.0051443
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course -0.0005526  0.0036129
##                                                        t value Pr(>|t|)    
## (Intercept)                                            -68.226   <2e-16 ***
## SeasonDay                                               -1.092    0.275    
## LocationUSA: Union Co, Bucknell Farm                    -0.709    0.478    
## LocationUSA: Union Co, Bucknell Ropes Course             0.196    0.844    
## SeasonDay:LocationUSA: Union Co, Bucknell Farm           1.028    0.304    
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course  -0.153    0.878    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2722 on 832 degrees of freedom
## Multiple R-squared:  0.01091,    Adjusted R-squared:  0.00497 
## F-statistic: 1.836 on 5 and 832 DF,  p-value: 0.1033
  • ANOVA(model) results:
#anova(Mass) #interaction term NS, SeasonDay p = 0.008244

#display results as nice table with kable
kbl(anova(Mass), digits = 2)  %>%
  kable_classic(full_width = F, html_font = "Cambria")
Df Sum Sq Mean Sq F value Pr(>F)
SeasonDay 1 0.49 0.49 6.66 0.01
Location 2 0.06 0.03 0.39 0.67
SeasonDay:Location 2 0.13 0.06 0.87 0.42
Residuals 832 61.66 0.07 NA NA
  • SeasonDay is significant (p = 0.01)

Conclusion

  • Use log-transformed mass in downstream analysis

  • Mass decreases over the season (potentially reflects using up energy reserves, desiccation).

  • Plot for multi-panel figure

p_body_mass_by_season_day <- ggplot(Adult_Final, aes(x=as.numeric(SeasonDay), y=log(Mass))) +
  geom_point(aes(color = Location), size = 1, alpha = 0.7) +
  geom_smooth(method = "lm", se=FALSE, color = "gray23") +
  theme_classic() +
  ylab("log(Body mass) (g)") +
  xlab("Season Day") +
  theme(plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90), legend.position = "none") +
  scale_color_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759")) + stat_regline_equation(label.x = 23, label.y = -2.8, aes(label = ..eq.label..)) +
  stat_regline_equation(label.x = 23, label.y = -2.6, aes(label = ..rr.label..))

p_body_mass_by_season_day

Length & Mass

Q1. Accounting for allometry (log-log), does elytron length predict body mass?

  • Summary
Both<-lm(log(Mass)~log(ElytralLength), data=Adult_Final)
summary(Both)
## 
## Call:
## lm(formula = log(Mass) ~ log(ElytralLength), data = Adult_Final)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.86176 -0.13766  0.00558  0.13144  0.66330 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        -7.28523    0.13514  -53.91   <2e-16 ***
## log(ElytralLength)  1.73265    0.06159   28.13   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1957 on 836 degrees of freedom
## Multiple R-squared:  0.4863, Adjusted R-squared:  0.4857 
## F-statistic: 791.4 on 1 and 836 DF,  p-value: < 2.2e-16
  • ANOVA
kbl(anova(Both), row.names = TRUE, digits = 130) %>%
  kable_classic(full_width = F, html_font = "Cambria")
Df Sum Sq Mean Sq F value Pr(>F)
log(ElytralLength) 1 30.31656 30.31656157 791.4157 4.705031e-123
Residuals 836 32.02444 0.03830675 NA NA
  • Plot
#using ggpubr package to get equations on graphs
#SOURCE: https://www.roelpeters.be/how-to-add-a-regression-equation-and-r-squared-in-ggplot2/
p_body_mass_by_elytron_length <- ggplot(Adult_Final,aes(x = log(ElytralLength), y = log(Mass))) + 
  geom_point(data = Adult_Final, aes(color = Location), size = 1, alpha = 0.7) + 
  geom_smooth(method = "lm", se=FALSE, color = "gray23") +
  theme_classic() + 
  stat_regline_equation(label.y = -2.75, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = -2.87, aes(label = ..rr.label..)) +
  scale_color_manual(name = "Location", labels = c("BNA", "BUF", "FDBCC"), values = c("#4E79A7", "#F28E2B", "#E15759")) +
  theme(legend.position = "none") +
  ylab("log(Body mass) (g)") +
  xlab("log(Elytron length) (mm)")

p_body_mass_by_elytron_length

Q2. If account for log(elytron length) first in the model, do Season Day and/or Location predict body mass?

  • Summary of effects:
#including interaction, but restricting to 2nd order interactions due to small sample size
Mass<-lm(log(Mass)~log(ElytralLength)*SeasonDay*Location, data=Adult_Final)
summary(Mass)
## 
## Call:
## lm(formula = log(Mass) ~ log(ElytralLength) * SeasonDay * Location, 
##     data = Adult_Final)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.81332 -0.13162  0.00501  0.12579  0.65748 
## 
## Coefficients:
##                                                                            Estimate
## (Intercept)                                                               -7.231751
## log(ElytralLength)                                                         1.698810
## SeasonDay                                                                 -0.032489
## LocationUSA: Union Co, Bucknell Farm                                       1.997806
## LocationUSA: Union Co, Bucknell Ropes Course                              -1.037070
## log(ElytralLength):SeasonDay                                               0.015145
## log(ElytralLength):LocationUSA: Union Co, Bucknell Farm                   -0.866447
## log(ElytralLength):LocationUSA: Union Co, Bucknell Ropes Course            0.520282
## SeasonDay:LocationUSA: Union Co, Bucknell Farm                            -0.013209
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course                     0.073983
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Farm          0.004831
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course -0.035764
##                                                                           Std. Error
## (Intercept)                                                                 0.848913
## log(ElytralLength)                                                          0.379069
## SeasonDay                                                                   0.056310
## LocationUSA: Union Co, Bucknell Farm                                        1.350131
## LocationUSA: Union Co, Bucknell Ropes Course                                1.054078
## log(ElytralLength):SeasonDay                                                0.025214
## log(ElytralLength):LocationUSA: Union Co, Bucknell Farm                     0.612366
## log(ElytralLength):LocationUSA: Union Co, Bucknell Ropes Course             0.474623
## SeasonDay:LocationUSA: Union Co, Bucknell Farm                              0.078929
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course                      0.061974
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Farm           0.035582
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course   0.027855
##                                                                           t value
## (Intercept)                                                                -8.519
## log(ElytralLength)                                                          4.482
## SeasonDay                                                                  -0.577
## LocationUSA: Union Co, Bucknell Farm                                        1.480
## LocationUSA: Union Co, Bucknell Ropes Course                               -0.984
## log(ElytralLength):SeasonDay                                                0.601
## log(ElytralLength):LocationUSA: Union Co, Bucknell Farm                    -1.415
## log(ElytralLength):LocationUSA: Union Co, Bucknell Ropes Course             1.096
## SeasonDay:LocationUSA: Union Co, Bucknell Farm                             -0.167
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course                      1.194
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Farm           0.136
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course  -1.284
##                                                                           Pr(>|t|)
## (Intercept)                                                                < 2e-16
## log(ElytralLength)                                                        8.46e-06
## SeasonDay                                                                    0.564
## LocationUSA: Union Co, Bucknell Farm                                         0.139
## LocationUSA: Union Co, Bucknell Ropes Course                                 0.325
## log(ElytralLength):SeasonDay                                                 0.548
## log(ElytralLength):LocationUSA: Union Co, Bucknell Farm                      0.157
## log(ElytralLength):LocationUSA: Union Co, Bucknell Ropes Course              0.273
## SeasonDay:LocationUSA: Union Co, Bucknell Farm                               0.867
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course                       0.233
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Farm            0.892
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course    0.200
##                                                                              
## (Intercept)                                                               ***
## log(ElytralLength)                                                        ***
## SeasonDay                                                                    
## LocationUSA: Union Co, Bucknell Farm                                         
## LocationUSA: Union Co, Bucknell Ropes Course                                 
## log(ElytralLength):SeasonDay                                                 
## log(ElytralLength):LocationUSA: Union Co, Bucknell Farm                      
## log(ElytralLength):LocationUSA: Union Co, Bucknell Ropes Course              
## SeasonDay:LocationUSA: Union Co, Bucknell Farm                               
## SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course                       
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Farm            
## log(ElytralLength):SeasonDay:LocationUSA: Union Co, Bucknell Ropes Course    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1933 on 826 degrees of freedom
## Multiple R-squared:  0.5048, Adjusted R-squared:  0.4982 
## F-statistic: 76.53 on 11 and 826 DF,  p-value: < 2.2e-16
  • ANOVA:
kbl(anova(Mass), row.names = TRUE, digits = 130) %>%
  kable_classic(full_width = F, html_font = "Cambria")
Df Sum Sq Mean Sq F value Pr(>F)
log(ElytralLength) 1 30.31656157 30.31656157 811.0826030 7.886660e-125
SeasonDay 1 0.16518932 0.16518932 4.4194385 3.583390e-02
Location 2 0.29247175 0.14623587 3.9123623 2.036480e-02
log(ElytralLength):SeasonDay 1 0.02514312 0.02514312 0.6726735 4.123577e-01
log(ElytralLength):Location 2 0.39495934 0.19747967 5.2833276 5.248466e-03
SeasonDay:Location 2 0.15246708 0.07623354 2.0395354 1.307438e-01
log(ElytralLength):SeasonDay:Location 2 0.12006939 0.06003470 1.6061550 2.012837e-01
Residuals 826 30.87414249 0.03737790 NA NA

Conclusions

  • log(Elytron length) sig. correlated with body mass, but can’t explain everything (v small p, r^2 = 0.49)

  • Location and SeasonDay significant when elytra length accounted for

  • Question: take residuals as “adult reserves” measurement?

    • no -> see Darlington, R. B., and T. V. Smulders. 2001. Problems with residual analysis. Anim. Behav. 62:599–602.

Combined multi-plot visualizing main results

#using patchwork
#get the plots together
multiplot <- p_elytral_length_by_season_day_BNA + p_elytral_length_by_season_day_BUF + p_elytral_length_by_season_day_FDBCC + p_body_mass_by_season_day + p_body_mass_by_elytron_length

#add annotations and layouts
w <- multiplot + plot_annotation(tag_levels = 'A') + plot_layout(widths = unit(c(3),c("in")), heights = unit(c(2), c("in")), guides = "collect") & theme(legend.position = 'bottom')

w

#ggsave("FigS1.EL_Mass_analysis.png", w, width = 12, height = 8, device = "png")

Step 4: Session Info

sessionInfo()
## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] patchwork_1.1.2  ggpubr_0.4.0     cowplot_1.1.1    sjPlot_2.8.11   
##  [5] kableExtra_1.3.4 gridExtra_2.3    fabricatr_1.0.0  ggthemes_4.2.4  
##  [9] qqplotr_0.0.5    Hmisc_4.7-1      Formula_1.2-4    survival_3.4-0  
## [13] lattice_0.20-45  dplyr_1.0.10     tidyr_1.2.1      ggplot2_3.3.6   
## 
## loaded via a namespace (and not attached):
##   [1] minqa_1.2.4         colorspace_2.0-3    ggsignif_0.6.3     
##   [4] deldir_1.0-6        ellipsis_0.3.2      sjlabelled_1.2.0   
##   [7] estimability_1.4.1  htmlTable_2.4.1     parameters_0.18.2  
##  [10] base64enc_0.1-3     rstudioapi_0.14     farver_2.1.1       
##  [13] fansi_1.0.3         mvtnorm_1.1-3       xml2_1.3.3         
##  [16] splines_4.2.1       cachem_1.0.6        robustbase_0.95-0  
##  [19] knitr_1.40          sjmisc_2.8.9        polynom_1.4-1      
##  [22] jsonlite_1.8.0      nloptr_2.0.3        ggeffects_1.1.3    
##  [25] broom_1.0.1         cluster_2.1.4       png_0.1-7          
##  [28] effectsize_0.7.0.5  compiler_4.2.1      httr_1.4.4         
##  [31] sjstats_0.18.1      emmeans_1.8.1-1     backports_1.4.1    
##  [34] Matrix_1.5-1        fastmap_1.1.0       cli_3.4.0          
##  [37] htmltools_0.5.3     tools_4.2.1         gtable_0.3.1       
##  [40] glue_1.6.2          Rcpp_1.0.9          carData_3.0-5      
##  [43] jquerylib_0.1.4     vctrs_0.4.1         svglite_2.1.0      
##  [46] nlme_3.1-159        insight_0.18.4      xfun_0.33          
##  [49] stringr_1.4.1       ps_1.7.1            lme4_1.1-30        
##  [52] rvest_1.0.3         lifecycle_1.0.2     rstatix_0.7.0      
##  [55] DEoptimR_1.0-11     MASS_7.3-58.1       scales_1.2.1       
##  [58] RColorBrewer_1.1-3  yaml_2.3.5          sass_0.4.2         
##  [61] rpart_4.1.16        latticeExtra_0.6-30 stringi_1.7.8      
##  [64] highr_0.9           bayestestR_0.13.0   checkmate_2.1.0    
##  [67] boot_1.3-28         rlang_1.0.5         pkgconfig_2.0.3    
##  [70] systemfonts_1.0.4   evaluate_0.16       purrr_0.3.4        
##  [73] htmlwidgets_1.5.4   labeling_0.4.2      tidyselect_1.1.2   
##  [76] processx_3.7.0      magrittr_2.0.3      R6_2.5.1           
##  [79] magick_2.7.3        generics_0.1.3      pillar_1.8.1       
##  [82] foreign_0.8-82      withr_2.5.0         mgcv_1.8-40        
##  [85] datawizard_0.6.0    abind_1.4-5         nnet_7.3-17        
##  [88] tibble_3.1.8        performance_0.9.2   modelr_0.1.9       
##  [91] car_3.1-0           interp_1.1-3        utf8_1.2.2         
##  [94] rmarkdown_2.16      jpeg_0.1-9          data.table_1.14.2  
##  [97] callr_3.7.2         digest_0.6.29       webshot_0.5.3      
## [100] xtable_1.8-4        munsell_0.5.0       viridisLite_0.4.1  
## [103] bslib_0.4.0